
% velocity of product pairs
% speed(10,9) = [19.075, 8.770]; %m/s, (10,9) pair
% speed(9,10) = [22.329, 10.266]; %m/s, (9,10) pair
% speed(10,10) = [16.559, 7.613]; %m/s, (10,10) pair
% speed(9,9) = [24.254, 11.151]; %m/s, (9,9) pair

% degeneracy of product pairs
% deg(10,9) = 29; % (10,9) pair 
% deg(9,10) = 29;  % (9,10) pair 
% deg(10,10) = 31;  % (10,10) pair 
% deg(9,9) = 28;  % (9,9) pair 
% speed(10,11)
%% simulation
    
c_10_9 = MC_coincidence(speed(10,9));

c_9_10 = MC_coincidence(speed(9,10));

c_9_9 = MC_coincidence(speed(9,9));

c_10_10 = MC_coincidence(speed(10,10));

ratios= [c_10_9 c_9_10 c_9_9 c_10_10].*[29 29 28 31]/sum([c_10_9 c_9_10 c_9_9 c_10_10].*[29 29 28 31])

   

function y = MC_coincidence(v)
    Omega = 20; 
    delta = 0;
    t = 0.02; %[us], relative delay between center of rising edge and 532 on 
    y = counts_rnd_angle(t,v,Omega,delta);
end


function ctot = counts_rnd_angle(mu,v,Rabi0,Delta0)  % measured counts at a certain pulse duration

    avg = 100000;

    cor = normrnd(0,1,[avg,3]);
    cor = cor./sqrt(sum(cor.^2,2));
    theta = acos(cor(:,3)); 
    phi = atan(cor(:,2)./cor(:,1));
    t0 = rand(1,avg)*50e-6; % with clean up, products only form in 50
%     r0 = v*t0;
    
    % assuming the dark mask inner radius is 125um
    t_dark = abs(125e-6./(v.*sin(theta))); %now in unit of s
    
    
    r0 = v.*t_dark+v.*t0';
    ctot=0;
    
    for i = 1:length(theta)
        th= theta(i);
        ph = phi(i);
        r0_vec = [r0(i,:)*sin(th)*cos(ph);r0(i,:)*sin(th)*sin(ph)];
        ctot = ctot+coi_counts(mu,v,th,ph,Rabi0,Delta0,r0_vec);
    end

end



function counts = coi_counts(mu,v,theta,phi,Rabi0,Delta0,r0)  % measured counts at a certain pulse duration
% assume each molecule only experience at most one REMPI pulse, which is
% valid for fast products

    rrabi_1 = Rabi_a(0,v(1)*sin(theta)*cos(phi),v(1)*sin(theta)*sin(phi),r0(:,1));
    rrabi_2 = Rabi_a(0,-v(2)*sin(theta)*cos(phi),-v(2)*sin(theta)*sin(phi),-r0(:,2));
    
    [p_e_1] = populations(rrabi_1,mu,v,theta,Rabi0,Delta0,1);  
    [p_e_2] = populations(rrabi_2,mu,v,theta,Rabi0,Delta0,2); 
   % p_c = sum(p_e_1.*p_e_2.*rrabi_1.^2.*rrabi_2.^2); % consider 532 power decay
    p_c = sum(p_e_1.*p_e_2); % ignore 532 power decay
    counts = p_c;
    
end



function [p_ee] = populations(Rabi,mu,v,theta,Rabi0,Delta0,ind) 
    f0 = [462906 445025]*1e3; % in unit of MHz
    Omega = Rabi0*Rabi;
    % detuning from Doppler shift
    Delta = v(ind)*cos(theta)/3e8*f0(ind)+Delta0;

    c0 = [Omega*2*pi 14*2*pi 0*2*pi Delta*2*pi 1 0 0];
    tlist = linspace(0,0.05,101); % [us]
    dt = 0.05/100;
    y0 = REMPI_Rate_solve(c0, tlist); % y(1) is rho00, y(2) is rho11, y(3) is rho10
    yy1 = y0(:,1); % ground state population
    yy2 = y0(:,2); % excited state population
    p_ee = yy2(1+round(mu/dt));
        
end
  


function R = Rabi_a(n,vx,vy,r0)
% r0 is vector (x0,y0)
% n is number of pulses 
% off_set is (center_r,theta), 
% no loss of generality to set theta as 0
    w = 0.75e-3; %[m] beam waist
    offset = 0;0.1e-3; 
    rx = r0(1)+n.*vx*100e-6 + offset;
    ry = r0(2)+n.*vy*100e-6;
    r = sqrt(rx.^2+ry.^2);
    if r>1e-3
        I = 0;
        R = 0;
    else
        I = exp(-2*r.^2./w.^2);
        R = exp(-1*r.^2./w.^2);
    end
end



function y_all = REMPI_Rate_solve(c, t)
    tspan = t;%[0 max(t)]; %[s]
    b = c(1:4);
    y0 = [c(5:7)];
    [~,y_all] = ode45(@(t,y) REMPI_RateEq(t,y,b), tspan, y0);
end

function dydt = REMPI_RateEq(t,y,b) % b are the coefficients
    % y(1) is rho00, y(2) is rho11, y(3) is rho10
    tmp = num2cell(b);
    [Omega01, Gamma1, GammaIon, Delta1] = tmp{:};

    dydt = zeros(3,1);
    dydt(1) = -1/2*1i*(Omega01*y(3)-conj(Omega01*y(3)));
    dydt(2) = -(Gamma1+GammaIon)*y(2)+ 1/2*1i *(Omega01*y(3)-conj(Omega01*y(3)));
    dydt(3) = -1/2*(Gamma1+GammaIon)*y(3)+ 1i*Delta1*y(3) + 1/2*1i *Omega01*(y(2)-y(1));
end


function ctot = deg(nK2, nRb2)
    Jtot = 1;
    Lmax = nK2 + nRb2 + 1;
    ctot = 0;
    for L = 0:Lmax
        for Np = abs(nK2 - nRb2):abs(nK2 + nRb2)
            for j4 = abs(Np - L):abs(Np + L)
                if mod(nK2 + nRb2 + L, 2) == 1 && j4 == Jtot
                    ctot = ctot + 1;
                end
            end
        end
    end
    % only care about allowed pairs
    E = 9.7711; % in unit of cm^-1
    BK2 = 5.478155*1e-2; % cm^-1
    DK2 = 7.8641*1e-8; % cm^-1
    BRb2 = 2.188943*1e-2; % cm^-1
    DRb2 = 1.29507*1e-8; %cm^-1
    U = BK2*nK2*(nK2+1)-DK2*(nK2*(nK2+1)).^2+BRb2*nRb2*(nRb2+1)-DRb2*(nRb2*(nRb2+1)).^2 ;
    T = E-U;
    if T<=0
        ctot=0;
    end
%     % manually imprint the 0.7 and 0.3 branching ratio between (e,o) and (o,e)
%     if mod(nK2,2)==0 && mod(nRb2,2)==1
%         ctot = ctot*0.7;
%     elseif mod(nK2,2)==1 && mod(nRb2,2)==0
%         ctot = ctot*0.3;
%     else 
%         ctot = 0;
%     end
        
end

function v = speed(nK2,nRb2)
    E = 9.7711; % in unit of cm^-1
    BK2 = 5.478155*1e-2; % cm^-1
    DK2 = 7.8641*1e-8; % cm^-1
    BRb2 = 2.188943*1e-2; % cm^-1
    DRb2 = 1.29507*1e-8; %cm^-1
    U = BK2*nK2*(nK2+1)-DK2*(nK2*(nK2+1)).^2+BRb2*nRb2*(nRb2+1)-DRb2*(nRb2*(nRb2+1)).^2 ;
    T = E-U;
    % constants
    h = 6.62606876*1e-34; %J*s
    c = 2.99792458*1e8; %m/s
    w2J = h*c*100;
    m0 =  1.66053906660*1e-27; % kg

    v_K2 = sqrt(2*87*2*T*w2J/80/(80+87*2)/m0);
    v_Rb2 = sqrt(2*40*2*T*w2J/(87*2)/(80+87*2)/m0);
    v = [v_K2,v_Rb2]; % in unit of m/s
end



